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(57) Abstract 

A parallel processing architecture and method, based on the 
3D Algebra Reconstruction Technique (ART), for iteratively recon- 
structing from cone beam projection data a 3D voxel data set (V) re- 
presenting an image of an object. The cone beam projection data is 
acquired by employing a cone beam x-ray source and a two-dimen- 
sional array detector to scan the object along a source scanning trajec- 
tory to obtain, at each of a plurality of discrete source positions (©j) 
on the source scanning trajectory, a 2D measured cone beam pixel 
data set The 3D voxel data set (V) is subdivided into or organ- 
ized as a plurality (m) of independent voxel subcubes (V°) through 
(V m ^) each containing a plurality of voxels. As a result of the sub- 
division of the 3D voxel data set (VHnto voxel subcubes, the 2D 
measured cone beam pixel data set 0^ (measured projection data 
array) is correspondingly subdivided for each source position (Qj) 
into projection data subsets, with overlapping regions between one 
or more adjacent projection data subsets corresponding to adja- 
cent voxel subcubes. Each voxel subcube and its corresponding 
projection data strip or subset is processed in parallel with other 
pairs of voxel subcubes and corresponding projection data strips 
or subsets, without interference. A bi-Ievel parallel-processor archi- 
tecture is employed for iterative reproject and merge, and split and 
backproject operations. 
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PARALLEL PROCESSING METHOD AND APPARATUS BASED ON THE ALGEBRA RBOONSTRUCTION 
TECHNIQUE FOR RECONSTRUCTING A THREE-DIMENSIONAL COMPUTERIZED TOMOGRAPHY 



Background of rh* Tnventton 

The present invention relates generally to three- 
dimensional (3D) computerized tomography (CT) and, more par- 
ticularly, to methods and apparatus employing parallel pro- 
cessing for reconstructing a 3D image of an object from cone 
1C beam projection data. 

In conventional computerized tomography for both 
medical and industrial applications, an x-ray fan beam and a 
linear array detector are employed. Two-dimensional (2D) 
imaging is achieved. While the data set is complete and 

15 image quality is correspondingly high, only a single slice of 
an object is imaged at a time. When a 3D image is required, 
a "stack of slices" approach is employed. Acquiring a 3D 
data set a 2D slice at a time is inherently tedious and time- 
consuming. Moreover, in medical applications, motion arti- 

20 facts occur because adjacent slices are not imaged simultane- 
ously. Also, dose utilization is less than optimal, because 
the distance between slices is typically less than the x-ray 
collimator aperture, resulting in double exposure to many 
parts of the body. 

25 A more recent approach, based on what is called 

cone beam geometry, employs a two-dimensional array detector 
instead of a linear array detector, and a cone beam x-ray 
source instead of a fan beam x-ray source, for much faster 
data acquisition. To acquire cone beam projection data, an 

30 object is scanned, preferably over a 360 # angular range, 

either by moving the x-ray source in a circular scanning tra- 
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jectory, for example, around the object, while keeping the 2D 
array detector fixed with reference to the source, or by 
rotating the object while the source and detector remain sta- 
tionary. In either case, it is relative movement between the 
5 source and object which effects scanning, compared to the 
conventional 2D "stack of slices" approach to achieve 3D 
imaging of both medical and industrial objects, with improved 
dose utilization. 

However, image reconstruction becomes complicated 
10 and requires massive computational power when a 3D image is 
reconstructed from cone beam projection data, in contrast to 
reconstruction of a 2D image from fan beam projection data. 

Literature discussing the cone beam geometry *for 3D 
imaging and generally relevant to the subject matter of the 

15 invention includes the following: Gerald N. Minerbo, "Convo- 
lutional Reconstruction from Cone-Beam Projection Data", IEEE 
Trans. Nucl. Sci., Vol. NS-26, No. 2, pp. 2682-2684 (April 
1979); Heang K. Tuy, "An Inversion Formula for Cone-Beam 
Reconstruction", SIAM J. Math., Vol. 43, No. 3, pp. 546-552 

20 (June 1983); L.A. Feldkamp, L.C. Davis, and J.W. Kress, 

"Practical Cone-Beam Algorithm", J. Opt. Soc. Am. A., Vol. 1, 
No. 6, pp. 612-619 (June 1984); Bruce D. Smith, "Image 
Reconstruction from Cone-Beam Projections: Necessary and 
Sufficient Conditions and Reconstruction Methods", IEEE 

25 Trans. Med. Imag., Vol. MI-44, pp. 14-25 (March 1985); and 

Hui Hu, Robert A. Kruger and Grant T. Gullberg, "Quantitative 
Cone-Beam Construction", SPIE Medical Imaging III: Image 
Processing, Vol. 1092, pp. 492-501 (1989). In general, this 
literature discloses various formulas for reconstruction of 

30 an image, including the use of a 3D Fourier transform or a 3D 
Radon transform. Questions of data completeness achieved 
with various source scanning trajectories are also con- 
sidered. 
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15 



20 



25 



'^^^El^^l^^^l^^^^^^^B^^^^^I The ART is an iter- 
ative algorithm which uses raysums to correct each voxel 
5 (volume element) of the reconstructed image at each itera- 
tion. The ART was initially disclosed in Richard Gordon, 
Robert Bender and Gabor T. Herman, "Algebraic Reconstruction 
Techniques (ART) for Three-dimensional Electron Microscopy 
and X-ray Photography", J. Theor. Biol., 29, pp. 471-481 
10 (1970) . Details regarding the mathematical foundations of 
the ART algorithm were described in Gabor T. Herman, Arnold 
Lent and Stuart W. Rowland, "ART: Mathematics and Applica- 
tions — A Report on the Mathematical Foundations and on the 
Applicability to Real Data of the Algebraic Reconstruction 
Techniques", J. Theor. Biol., 43, pp. 1-32 (1973). Extension 
of the additive ART to 3D cone beam geometry is described in 
M. Schlindwein, "Iterative Three-Dimensional Reconstruction 
from Twin-Cone Beam Projections", IEEE Trans. Nucl. Sci., 
Vol. NS-25, No. 5, pp. 1135-1143 (October 1978). The ART was 
recently used for vascular reconstruction, as reported in A. 
Rougea, K.M. Hanson and D. Saint -Felix, "Comparison of 3-D 
Tomographic Algorithms for Vascular Reconstruction", SPIE, 
Vol. 914 Medical Imaging II, pp. 397-405 (1988). 

In general, the Algebra Reconstruction Technique 
requires a large amount of computation time, and has not been 
implemented with parallel processing techniques. 

In view of the parallel processing aspect of the 
present invention, the system described in Richard A. Robb, 
Arnold H. Lent, Barry K. Gilbert, and Aloysius Chu, "The 
30 Dynamic Spatial Reconstructor", J. Med. Syst., Vol. 4, No. 2, 
pp. 253-288 (1980) is relevant. The Dynamic Spatial Recon- 
structor employs twenty-eight x-ray sources and twenty-eight 
x-ray imaging systems in a synchronous scanning system to 
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acquire data all at one time for a subsequent "stack of 
slices" reconstruction using conventional 2D reconstruction 
algorithms. The Robb et al. reference refers to the use of 
"high-speed parallel processing techniques" for the recon- 
5 struction computation. 

Summary of thf> Tnw>nf^ n 

Accordingly, it is an object of the invention to 
provide practical methods and apparatus for reconstructing 3D 
images from cone beam projection data, in view of the massive 
computational power required. 

10 It is another object of the invention to provide 

methods and apparatus employing parallel processing for 
reconstructing 3D images from cone beam projection data. 

It is another object of the invention to provide 
methods and apparatus which can reconstruct selected regions 
15 of interest/ without spending processor time reconstructing 
regions outside the regions of interest. 

Briefly, and in accordance with an overall aspect 
of the invention, a parallel processing architecture and 
method, based on the 3D Algebra Reconstruction Technique 

20 (ART) , are defined for iteratively reconstructing from cone 
beam projection data a 3D voxel data set V representing an 
image of an object. The cone beam projection data is 
acquired by employing a cone beam x-ray source and a two- 
dimensional array detector to scan the object along a source 

25 scanning trajectory to obtain, at each of a plurality of 
discrete source positions 6* on the source scanning tra- 
jectory, a 2D measured cone beam pixel data set P; . The 
invention is applicable to reconstructing an object image 
from cone beam projection data acquired with any scanning 

30 geometry, so long as it is well defined. Examples include 
source scanning trajectories which comprise a single circle, 
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dual or multiple parallel circles, dual perpendicular cir- 
cles, or square waves on a cylindrical surface. However, in 
the embodiment described in detail hereinbelow, a circular 
source scanning trajectory lying in a plane perpendicular to 
5 a rotation axis passing through the object and perpendicular 
to all planes of the detector is employed. In this case, the 
discrete source positions 9/ are defined simply as angular 

positions. With various other scanning trajectories, the 
discrete source positions 8/ are defined by parameters in 

10 addition to or instead of angular position. 

In accordance with the invention, the 3D voxel data 
set V is subdivided into or organized as a plurality of m of 
voxel subcubes V° through V**- 1 which are independent of each 
other. In other words, there is no overlap between adjacent 

15 voxel subcubes. Each voxel subcube contains a plurality of 
voxels. In the illustrated embodiment, boundaries between 
adjacent voxel subcubes lie in planes which are perpendicular 
to the detector planes and parallel to each other. However, 
independent voxel subcubes may be organized in a variety of 

20 ways. 

As a result of the subdivision of the 3D voxel data 
set V into voxel subcubes, the 2D measured cone beam pixel 
data set P t (measured projection data array) is correspond- 
ingly subdivided for each source position 8/. In the embodi- 

25 ment described in detail herein wherein boundaries between 
adjacent voxel subcubes lie in parallel planes, the measured 
projection data array is divided into what may be termed pro- 
jection data strips. In the more general case, the measured 
projection data array is divided into what may be termed pro- 

30 jection data subsets. 



Thus, for each source position 8/, each voxel sub- 
cube has a corresponding projection data strip or projection 
data subset. However, unlike the voxel subcubes, the projec- 
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tion data strips or subsets are not independent. In other 
words, dependent upon the particular scanning geometry 
defined and the particular voxel and pixel positions for each 
particular source position 0/, there are overlapping regions 

5 between one or more adjacent projection data strips or sub- 
sets corresponding to adjacent voxel subcubes. 

Nevertheless, in accordance with the invention, 
each voxel subcube and its corresponding projection data 
strip or subset can be processed in parallel with other pairs 

10 of voxel subcubes and corresponding projection data strips or 
subsets, without interference. A bi-level parallel-processor 
architecture is employed. Backpro ject/repro ject processor in 
level 0 selectively perform voxel driven backprojection from 
the projection data strips or subsets to the corresponding 

15 voxel subcubes, and voxel-driven reprojection from the voxel 
subcubes to the corresponding projection data strips or sub- 
sets. Each backpro ject/reproject processor operates on data 
which is independent of the data for the other backpro- 
ject/reproject processors, so the processors can operate in 

20 parallel. At least one split/merge processor in level 1 

splits a projection data set for a particular source position 
8/, into projection data strips or subsets from reprojection 

into a projection data set. During iterative reconstruction 
implementing the Algebra Reconstruction Technique, for each 
25 source position 6/, data are sent back and forth between level 

0 for backprojection and reprojection and level 1 for splitt- 
ing and merging. 

Preferably, level 1 is organized in a tree-type 
manner as a plurality of layers. The split/merge processors 
30 can comprise transputers. 

More particularly, a method in accordance with the 
invention includes the steps of organizing the voxel data set 
V as a plurality m of independent voxel subcubes V° through 
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V™- 7 , each voxel subcube including a plurality of voxels; ini- 
tializing the voxel data set V; and performing successive 
iterations to correct the values of at least selected voxels 
of the 3D voxel data set based on the measured cone beam data 

A 

5 sets P s . In a particular embodiment, boundaries between adja- 
cent voxel subcubes lie in planes perpendicular to all planes 
of the 2D detector array and parallel to each other. 

During each iteration/ the values of at least the 
selected voxels for each of the discrete source positions 6/ 

10 on the source scanning trajectory are corrected by reproject- 
ing each voxel subcube V° through V™' 1 to calculate pixel val- 
ues for a corresponding 2D calculated projection data strip 
or subset from a group of calculated projection data strips 
or subsets If through /J""" 1 of a 2D calculated projection data 

15 set P± including a plurality of pixels , the projection data 
strips of subsets through P"~ x overlapping in part; merging 
the calculated projection data strips of subsets P? through 
P*~ l to calculate pixel values of the 2D calculated projection 
data set Pi, the step of merging including adding pixel values 

20 in any overlapping regions of the calculated projection data 
strips of subsets tf through /J* -1 ; calculating a 2D correction 
projection data set £/ including a plurality of pixels by 
determining the normalized difference between the value of 
each pixel of the measured projection data set P and the 

25 value of the corresponding pixel of the calculated projection 
data set Pi splitting the 2D correction projection data set £,* 

into a plurality of 2D correction projection data strips or 
subsets £f through £/"~ l corresponding to the voxel subcubes V° 
through V**' 1 , the 2D correction projection data strips or sub- 
30 sets Ef through £*"* overlapping in part, including duplicat- 
ing element values for any overlapping regions of the correc- 
tion projection data strips or subsets £/ through £*" 1 ; and 
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backprojecting each correction projection data strip or sub- 
set £/ through E?~ x to correct voxel values in the correspond- 
ing voxel subcube of the plurality of voxel subcubes V° 
through V**' 1 . 

5 The step of merging is organized in a tree-type 

manner as a plurality of layers and includes, for each layer , 
merging groups of calculated projection data strips or sub- 
sets from the next lower layer to present a lesser number of 
calculated projection data subsets to the next higher layer, 
10 the lowest layer merging groups of the calculated projection 
data strips or subsets Ff through P"" 1 from the step of repro- 

jecting, and the highest layer producing the calculated pro- 
jection data set Pf. 



The step of splitting is also organized in a tree- 
15 type manner as a plurality of layers and includes, for each 
layer, splitting one or more correction projection data 
strips or subsets from the next higher layer to present a 
greater number of correction projection data strips or sub- 
sets to the next lower layer, the highest layer splitting the 
20 correction data set £/, and the lowest layer producing the 

correction projection data strips or subsets E° through £*" 1 . 



The method further includes employing a plurality m 
of backproject/reproject processors corresponding to the 
voxel subcubes to perform the steps of reprojecting and back- 
25 projecting, and employing a plurality of split/merge proces- 
sors connected in a tree structure to perform the steps of 
merging and splitting. 

The method additionally includes, during each iter- 
ation, for each of the discrete source positions 0/, in con- 

30 junction with the steps of reprojecting each voxel subcube V° 
through V™- 1 , calculating data values for a corresponding 2D 
normalizing factor data strip or subset from a group of 2D 
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normalizing factor strips or subsets N? through Af"" 1 of a 2D 
normalizing factor data set Aft including a plurality of pix- 
els, the 2D normalizing factor strips or subsets N? through 
overlapping in part; merging the normalizing factor data 
5 strips or subsets AT* through Aft*" 1 to calculate element values 
of the 2D normalizing factor data set Aft; and f during the step 
of calculating the 2D correction data set £/, employing the 
value of each corresponding element of the normalizing factor 
data set Aft. 

10 The step of reprojecting is voxel driven and 

includes, for each voxel subcube of the plurality of voxel 
subcubes V° through W 1 ' 1 , initializing the corresponding cal- 
culated projection data strip or subset from the group of 
projection data subsets P° through /J"* 1 ; initializing the cor- 

15 responding normalizing factor data strip or subset from the 
group of normalizing factor data strips or subsets N° through 
Aft"" 1 ; and for each of the selected voxels of the 3D voxel data 

set V included in the particular voxel subcube, adding to the 
value of each pixel of the corresponding calculated projec- 

20 tion data strip or subset from the group of calculated pro- 
jection data strips or subsets P? through influenced by 
the selected voxel the voxel value multiplied by a weighting 
factor h(p,v) determined for the voxel and pixel positions 
based on geometry such that pixel values are cumulated in the 

25 corresponding artificial projection data subset, and adding 
to the value of each element of the corresponding normalizing 
factor data strip or subset from the group of normalizing 
factor data strips or subsets N? through N?~ l the square of 
the weighting factor h(p,v) determined for the voxel and pixe 

30 positions based on geometry such that element values are 

cumulated in the corresponding normalizing factor data strip 
or subset. 
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The step of backprojecting is also voxel driven and 
includes , for each voxel subcube of the plurality of voxel 
subcubes V° through W* 1 , for each of the selected voxels of 
the 3D voxel data set V included in the particular voxel sub- 
5 cube, adding to the voxel value the value of each element of 
the corresponding 2D correction projection data strip or sub- 
set of the group of correction projection data strips or sub- 
sets £* through £*~ l in a pixel position influenced by the 
selected voxel multiplied by a weighting factor h(p,v) deter- 
10 mined for the voxel and pixel positions based on geometry 
such that corrected voxel values are cumulated* 

Parallel processing apparatus implementing the 
Algebra Reconstruction Technique in accordance with the 
invention includes a data memory for storing the voxel data 

15 set V organized as a plurality m of independent voxel sub- 
cubes V° through V»»-J, each voxel subcube including a plural- 
ity of voxels; and processing elements for initializing the 
voxel data set V and performing successive iterations to cor- 
rect the values of at least selected voxels of the 3D voxel 

20 data set based on the measured cone beam data sets P it during 

each iteration correcting the values of at least the selected 
voxels for each of the discrete source positions 8/ on the 
source scanning trajectory . 

The processing elements include in a level 0 a plu- 
25 rality m of backpro ject/repro jectors corresponding to the 
voxel subcubes, and in a level 1 at least one split /merge 
processor. The backpro ject/repro ject processors are operable 
to reproject each voxel subcube V° through V™' 1 to calculate 
pixel values for a corresponding 2D calculated projection 
30 data subset from a group of calculated projection data sub- 
sets P? through P~~ x of a 2D calculated projection data set Pi 
including a plurality of pixels, the projection data subsets 
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P? through P"~ x overlapping in part. The split/merge proces- 
sor is operable to merge the calculated projection data sub- 
sets P? through P"~ l to calculate pixel values of the 2D 
calculated projection data set Pi, adding pixel m values in any 
5 overlapping regions of the calculated projection data subsets 
ff through P*~ l . The processing elements include means for 
calculating a 2D correction projection data set £i including 
a plurality of pixels determining the normalized difference 
between the value of each pixel of the measured projection 
10 data set pixel P : and the value of the corresponding pixel of 
the calculated projection data set The split/merge pro- 

cessor is also operable to split the 2D correction projection 
data set £/ into a plurality of 2D correction projection data 

subsets £* through E?~ l corresponding to the voxel subcubes V° 
15 through V 971 ^, duplicating element values for any overlapping 
regions of the correction projection data subsets £f through 

£*~ l . The backpro ject/repro ject processors are also operable 
to backproject each correction projection data subset £/ 
through £*"* to correct voxel values in the corresponding 
20 voxel subcube of the plurality of voxel subcubes V° through 

Preferably there are a plurality of split/merge 
processors organized as a tree-type structure in a plurality 
of layers . The split/merge processors are operable to merge 
25 groups of calculated projection data subsets from the next 
lower layer to present a lesser number of calculated projec- 
tion data subsets to the next higher layer, the lowest layer 
merging groups of the calculated projection data subsets Pf 
through P"~ x from the backpro ject/repro ject processors , and 

30 the highest layer producing the calculated projection data 

set Pi. The split /merge processors are also operable to split 

one or more correction projection data subsets from the next 
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higher layer to present a greater number of correction 
projection data subsets to the next lower layer, the highest 
layer splitting the correction data set £/, and the lowest 

layer producing the correction projection data subsets £° 

5 through £/*~ l . 

At least the split/merge processors comprise trans- 
puters. The split/merge processors of each layer each 
include on I/O connected to the next higher layer, and three 
I/O parts connected to the next lower layer. 

Brifrf Description of th* Draw^gs 

While the novel features of the invention are set 
forth with particularity in the appended claims, the inven- 
tion, both as to organization and content, will be better 
understood and appreciated, along with other objects and fea- 
tures thereof, from the following detailed description taken 
in conjunction with the drawings, in which: 

FIG. 1 depicts data acquisition employing cone beam 
geometry for scanning an object to acquire cone beam projec- 
tion data, and additionally graphically depicts the geometry 
for backprojection and reprojection operations; 

20 FIG. 2 represents the steps in the prior art 3D 

Algebra Reconstruction Technique (ART) for reconstruction; 

FIG. 3 depicts the organization of a 3D voxel data 
set into independent voxel subcubes; 

FIG* 4 depicts the geometry relationship between 
25 voxel subcubes and projection strips or subsets; 

FIG. 5 is a black diagram of apparatus implementing 
ART in accordance with the invention; 
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FIG. 6 represents the steps of reprojection and 
merge in the ART implementation of the invention; 

FIG. 7 represents the steps of split and backpro- 
ject in the ART implementation of the invention; 

5 FIG. 8 depicts the merging of projection strips or 

subsets; and 

FIG. 9 depicts splitting into projection strips or 

subsets. 

Detailed Description 

Referring first to FIG. 1, depicted in a typical 
10 scanning and data acquisition configuration employing cone 
beam geometry. A cube 20 alternatively represents an actual 
3D object being scanned or a 3D voxel (volume element) data 
set V representing the actual object. The voxel data set V 
may be indexed, for example, by voxel number v. It will be 
15 appreciated that, during a scanning and data acquisition 
operation, the actual object is scanned with x-rays; and, 
during a subsequent image reconstruction operation, the voxel 
data set V is calculated to represent an image of the actual 
object. The present invention is directed to image recon- 
20 struction. 

A global coordinate system C] = (X ¥f Y v9 Z v ) is defined, 
with its origin (0,0,0) at the center of the voxel data set V 
or object 20. The object 20 is positioned within a field of 
view between a cone beam x-ray point source 22 and a 2D array 
25 detector 24, from which cone beam projection data is 

acquired. The 2D array detector 24 is connected to a conven- 
tional high speed data acquisition system (not shown). An 
axis of rotation 26 coincident with the Z v axis passes through 

the field of view and the object 20. Perpendicular to the Z v 
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axis of rotation 26 is a midplane within which the X v and Y v 

axes lie, as well as a circular source scanning trajectory 28 
centered on the Z v axis of rotation 26. 

For scanning the object at a plurality of discrete 
5 source position 8i on the source scanning trajectory 28, 

which in FIG. 1 are discrete angular positions 0/, the x-ray 
source 22 moves relative to the object 20 along the source 
scanning trajectory 28, while the 2D array detector 24 
remains fixed in position relative to the source 22. Either 
10 the source 22 and detector 24 may be rotated around a sta- 
tionary object 20, or the object 20 may be rotated while the 
source 22 and detector 24 remain stationary. 

At each discrete source position 9/ a 2D measured 
cone beam pixel data set P g is acquired by the data acquisi- 
15 tion system (not shown) connected to the array detector 24 , 
and stored for subsequent reconstruction. Each measured cone 
beam pixel data set P t includes a plurality of pixels indexed, 
for example, by pixel number p. 



20 



Thus, in the particular geometry of FIG. 1, data 
are acquired at a number of angular positions or view angles 
Qi around the object. As depicted, a is the source-to-detec- 
tor distance and j$ is the source-to-rotation-axis distance. 
For each view angle 8i, the x-ray source 22 is located at 

[psin9 i9 -ficosd i$ 0Y and a 2D projection P* is acquired. The cen- 
25 ter of the 2D projection array 24 is located at 
0, =[-(a-j3)sin0 |f (a-/J)cose,,of . 

FIG. 2 summarizes the procedure of the prior art 
three-dimensional Algebra Reconstruction Technique (ART) , 
which is an iterative procedure using raysums to correct the 
30 value of each voxel of the 3D voxel data set V at each itera- 
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tion. After a number of iterations, the values of the voxel 
data set.V converge to represent the actual 3D object. 

In box 30, the procedure commences by initializing 
the 3D voxel data set V, for example by clearirig all voxel 
5 values to zero. The procedure then goes through successive 
iterations and, for each iteration, through each source posi- 
tion or projection angle Q\ corresponding to discrete angular 
positions on the source scanning trajectory 28. For each 
projection angle 0/ there are three steps: repro jection, 

10 error calculation and backpro jection. The procedure is com- 
pleted (not shown) either when a predetermined number of 
iterations have been processed, or when the error is below a 
predetermined threshold. The projection angles 8/ may be 

processed in any order, either sequentially or in some other 
15 order Depending upon the particular scanning geometry, the 

order of projection angle 0/ processing can affect the rate of 

convergence of the voxel data set V. 

Box 32 states the step of repro jection, which is in 
turn detailed in Box 34 in the manner of a subroutine. The 

20 3D voxel data set V (or a selected portion of interest) is 

reprojected from real space to projection space to generate a 
2D calculated projection data set Pi for the particular pro- 
jection angle 8/. The 2D calculated projection data set P± 
corresponds on a pixel-by-pixel basis with the actual mea- 

25 sured cone beam pixel data set P., but with actual pixel data 
values differing depending upon how accurately the voxel data 
set V represents the actual object, and depending upon vari- 
ous inherent computational errors. In conjunction with the 
reprojection, a 2D normalizing factor data set Ni is also 

30 calculated, corresponding on a pixel-by-pixel basis with the 
calculated and measured projection data sets Pi and P : . 

Considering the reprojection process of Box 34 in 
detail, the calculated projection data set Pi and the normal- 
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izing factor data set Nf are each initialized by clearing all 

pixels to zero. Then for each voxel v of the 3D voxel data 
set V, each pixel p of the calculated projection data set P/ 
and the corresponding normalizing factor data set Ni influ- 

5 enced by the voxel is identified in a manner described here- 
inbelow with reference again to FIG. 1. For each pixel p so 
influenced by voxel v, the following three operations are 
performed: (1) A weighting factor h(p,v) for the particular 
voxel and pixel positions is determined based on geometry, 
10 for example by a well-known bi-linear interpolation process. 
(2) The value of the voxel V(v) is multiplied by the weight- 
ing factor h(p,v) and added to the value of the calculated 
projection data set pixel Pj(p) such that pixel values are 

cumulated. (3) The weighting factor h(p,v) is multiplied by 
15 itself and added to the value of the normalizing factor dat 
set element Nf(p) such that element values are cumulated. 

Returning to Box 30 of FIG. 2, a 2D correction pro- 
jection data set £i indexed, for example, by pixel number p, 
is calculated. For each pixel of the corresponding data 
20 sets, the difference between the value of the measured pro- 
jection data set pixel P^p) and the value of the corresponding 
calculated projection data set pixel Pi(p) is determined, and 

then normalized by dividing by the value of the corresponding 
normalizing factor data set element Ntfp). The result is the 

25 value of the corresponding correction projection data set 
element Etfp). 

Box 36 states the step of backprojection, which is 
in turn detailed in Box 38 in the manner of a subroutine. 
The 2D correction projection data set £/ for the particular 
30 projection angle 8/ is backpro jected from projection space to 

real space to correct voxel values in the 3D voxel data set 
V. 
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Considering the backprojection process of Box 38 in 
detail, for each voxel v of the 3D voxel data set V, each 
pixel p of the correction projection data set £,• influenced by 

the voxel is identified. For each pixel p so influenced by 
5 voxel v, the weighting factor h(p,v) is determined via inter- 
polation based on geometry, and then the value of the correc- 
tion projection data set pixel Etfp) is multiplied by the 
weighting factor h(p,v), and the product added to the value 
of the voxel value V(v), such that voxel values are cumula- 
10 tively corrected. 

Referring again to FIG. 1, the manner in which 
array detector pixels influenced by each voxel is determined 
will now be described in the context of the backprojection 
and reprojection procedures. 

Let V(x v , y v , z v ) represent a three-dimensional voxel 
data set defined in the three-dimensional discrete Cartesian 
coordinate system, C\, with the origin located at the center 
of the 3D cubic data. For a 3D voxel data set with 
(N%xN%xNZ) voxels, every voxel can be indexed by (/„/,,/,). 
The coordinate of each voxel is stated as 

y.=i,xK, / >= — x+i t ...,-JL_i (1) 

N v N v 
Z,=i t XK t / f =-i^ + l,...,i^-l 

where K x , K y , and K 2 denote the voxel resolution along the x, 
y, and z axis in the coordinate system C\, respectively. A 

two-dimensional discrete Cartesian coordinate system is 

25 used to address each detector element in the 2D area detec- 
tor. Every detector element DUd, yd) is indexed by 
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(2) 



the Ji. x and Jl y are used to indicate the detector resolution 
along the x and y axis, respectively. The relative relation- 
ship between and C£ can be stated as 4x4 transformation 
function T p ^ v . For each detector P $s (x d ,y d ) defined in the 
detector coordinate system C^, the coordinate <x d , y d ) can be 

transferred to the voxel coordinate system ^ by applying the 
transformation function ^ , i.e., 



(3) 



V 




V 


y 


-[v.] 


y* 






,0, 



10 



transformation function T F , which transfers the detector 

coordinate system to the voxel coordinate system, is 
expressed as rotating about the y axis by 8/, then rotating 

90* around the reference x axis, and finally translating by 



i.e. 
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(4) 



where 0/ is View angle measured respect to the y axis in C\ 

a is Source to detector distance (STOD) 

P is Source to rotation axis distance (STOD) 

then Equation (3) can be rewritten as 

'cos0 0 sin^-Ca-Zftsinfli 



sin0,. 0 -cos0,(a-/J)cos0; 
0 10 0 
0 0 0 



(5) 



i.e. 



x - x d cos 0, - (a - /J) sin 0; 
y = x d sin 6 { + (a - /?) cos 0 f 

. z =y<i 



(6) 



For each voxel V(x v ,y v ,z v ), given a view angle 9/, an integra- 
10 tion line such as an exemplary integration line 44 passes 
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through the x-ray source 22 located at [Psin& n -/3 cos0„O] r and 
voxel V(x v ,y v ,z v ) , can be expressed as 

x-PsmS, = y+Pcosd : = z • (?) 

where (x,y,z) is any point belonging to the integration line 

5 44. This integration line 44 intersects the detector plane 
at the location (x d ,y d ) in the projection coordinate system 

C£. . Substituting Equation (6) into Equation (7) yields 

^cosfl i -(g-/?)sing f -fisinfl £ _ x rf sing f + (a-j8)cosg g +ffcosg f 

x v -0siii0 2 y,+/Jcos0,. " z ¥ 

Then the intersection point (xd/ yd) on the detector or pro- 
10 jection plane can be calculated as 

x - ^cosB^y^mO,] 
4 ~ y.cosfl,. -Jc v sin0 f +0 



y v cos^-jr v sin^+j3 

A voxel-driven approach is applied for both back- 
projection and repro jection. At a given view angle 8,-, a 

15 unique intersection point <xd, yd) on the detector plane 

can be calculated for a given voxel (x v , y v , z v ) in C) , as 

represented in FIG. 1 by the intersection point 40 and voxel 
42 along line 44. Four immediate neighbor pixels surrounded 
the intersection point (xdr yd) • By a known process of bi- 
20 linear interpolation , weighting factors reflecting the con- 
tribution of each voxel to each pixel on the detector plane 
can be calculated for reprojection, and weighting factors 
reflecting the contribution of each pixel on the detector 
plane to each voxel can be calculated for backpro jection . 
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With reference now to FIG. 3, the 3D voxel data set 
V (or 20 in FIG. 1) is organized as a plurality m of indepen- 
dent voxel subcubes V° through V m "- ? . The voxel subcubes are 
independent in the sense that no overlapping region exists 
5 between adjacent voxel subcubes. While the voxel data set 
may be divided in a variety of ways, preferably boundaries 
between adjacent voxel subcubes lie in planes which are per- 
pendicular to the detector planes and parallel to each other, 
as depicted in FIG. 3. The size of each voxel subcube is 
10 constrained by the memory space associated with each proces- 
sor element used to do back projection and repro jection, 
described hereinbelow with reference to FIGS. 5, 8 and 9. 

FIG. 4 depicts the geometric relationship between 
the voxel subcubes of FIG. 3 and projections on the detector 
15 array for each angular position 6/. For each voxel subcube in 

real space there is a corresponding projection subset on the 
projection plane. In the planar arrangement of voxel sub- 
cubes illustrated, the projection subsets take the form of 
parallel strips on the detector plane. 

Thus, in FIG. 4 there are two representative voxel 
subcubes 4 6 and 48, with corresponding projection strips 50 
and 52. There is an overlapping region 54 between the two 
projection strips 50 and 52. It is however a feature of the 
invention that each voxel subcube and its corresponding pro- 
jection strip can be projected and reprojected in parallel 
with other pairs of voxel subcubes and projection strips 
without interference. 

FIG. 5 depicts apparatus in accordance with the 
invention, comprising a bi-level parallel processor architec- 
30 ture. Included in the FIG. 5 apparatus is a data memory 60 
for storing the 3D voxel data set V, organized as indicated 
into separate memory areas for each of the independent voxel 
subcubes V° through V* 1 ** . 



20 



25 



WO 92/05507 



-22- 



PCT/US91/06822 



In a level 0 are a plurality of m of backpro- 
ject/reproject processors, such as processor 62, correspond- 
ing respectively to the voxel subcubes . Each of the backpro- 
ject/reproject processors of level 0 operates independently 
5 of an in parallel with the other processors of level 0 to 

reproject and backproject back and forth between a particular 
voxel subcube and a particular projection strip as part of 
the overall iterative process. 

In a level 1 is at least one, and preferably many, 
10 split/merge processor, such as processor 64. Using 4-link 
connectable processors, such as the INMOS T800 transputer, 
processors in level 1 are organized in layers in a tree-type 
manner. Each level 1 processor (or node) has three sons and 
one parent in the particular embodiment illustrated. The 
15 INMOS T800 transputer is basically a computer on a chip with 
a 32-bit integer processor and a 64 bit floating point pro- 
cessor, with a practically achievable floating point computa- 
tion rate of approximately 1.5 MFLOPS. It includes 4 kBytes 
of on-chip memory, and address space for 4 Gbytes of off-chip 
20 memory. Each transputer can be connected to four neighbors 
over a point-to-point link with a 20-Mbit /second bandwidth. 
An available development system, namely a Meiko Computing 
Surface, includes 512 transputers, with a maximum computation 
rate of approximately 800 MFLOPS. 

25 The level 1 processors are operable, following a 

level 0 reproject operation, which produces m calculated pro- 
jection data subsets Ff through P"" 1 , to merge the calculated 

projection data subsets /f through P"~ l into a memory 66, as 
well as to merge normalizing factor data subsets Nf through 
30 iV/"~ l from the level 0 processors into a single normalizing 

factor data set Aft for storing in the memory 66. The level 1 
processors are also operable to split a correction projection 
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data set £/ stored in the memory 66 into a plurality of cor- 
rection projection data subsets £? through E?~ x for backpro- 

jection by the level 0 processors. In FIG. 5, the arrows 68 
and 70 respectively indicate these two types of operation, 
5 namely, reproject then merge, and split then backpro ject. 

Also shown in FIG. 5 is a correction processor 72 
which operates in a conventional manner as described herein- 
above to calculate the correction projection data set £ ( * based 

on the difference between the calculated projection data set 
10 Pi, and the measured projection data set P s for the particular 

projection angle 9i. The measured projection data sets 
resulting from data acquisition as described hereinabove with 
reference to FIG. 1 are stored in exemplary memories 74, 76, 
78 and 80. 

15 Thus,, one projection array is associated with each 

split/merge processor in level 1. During backpro jection, 
each split/merge processor reads in a projection reads in a 
projection strip from its parent node, splits it into three 
projection substrips, and then distributes those to its son 

20 nodes for further splitting or backpro jection . During repro- 
jection, each split/merge processor reads in three projection 
strips from its son nodes (either the backpro ject/repro ject 
processors or split/merger processors) , merges them into one 
projection strip, and sends that projection strip to its par- 

25 ent node for further merging. The number of processor layers 
in level 1 is related to the total number of backpro- 
ject/reproject processors in level 0. One voxel subcube and 
its corresponding projection strip are associated with each 
backpro ject/repro ject processor in level 0. During backpro- 

30 jection, each split/merge processor reads in a projection 
strip from its parent node and splits it into three projec- 
tion substrips, and then distributes those to its son nodes. 
Each backproject/reproject processor receives a projection 
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strip from its parent node and performs voxel-driven backpro- 
jection. During repro jection, each backpro ject/reproject 
processor calculates a projection strip for the desired view 
angle and sends the resulting projection strip to its parent 
5 node (split/merge proces^r) for merging. The 3D algebra 
Reconstruction Technique is achieved by iteratively perform- 
ing backpro jection and repro jection. The entire parallel pro- 
cessing procedure of the 3D cone beam algebra reconstruction 
algorithm can be divided into four major tasks: split/ 
10 merge, backpro ject, and reproject. 

These four tasks are represented in Boxes 82 and 84 
of FIGS . 6 and 7 . In the ART implementation of the inven- 
tion, the repro jection and merge tasks in Box 82 of FIG. 6 
are substituted for the prior art repro jection step repre- 
15 sented in Box 34 of FIG. 2. Similarly, the split and back- 
project tasks in Box 84 of FIG. 7 are substituted for the 
prior art backpro jection step represented in Box 38 of FIG. 
2. The four tasks will now be considered in detail in the 
order presented in the Boxes 82 and 84. 

20 The reproject process is executed by each backpro- 

ject/reproject processor in level 0 or repro jection. Its 
primary function is to generate a projection strip for a 
given view angle from a voxel subcube and send the results to 
its parent node for projection merging. The repro jection 

25 procedure uses a voxel-driven approach and is illustrated as 
follows : 
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for a given view angle 8j 

set P - 0 and N = 0 

for each voxel v in the voxel subcube 

calculate pixel location in the projection strip influenced by v for each 
5 pixel p influenced by v 

calculate h(v,p), the weighting factor, via interpolation 
P(p) = P(p) + h(p,v)V(v) 
N(p)«N(p) + h(p,v)h(p,v) 

end for 

10 end for 

end for 

The merge process runs on the split/merge proces- 
sors in level 1 for reprojection purposes. Its primary func- 
tion is to merge three projection strips, read from its three 
son nodes, into one. As mentioned hereinabove, an overlap- 
ping region between two projection strips corresponds to two 
adjacent voxel subcubes. The merge is achieved by adding 
projection strips together in the overlapping region accord- 
ing to its global coordinate system, as represented in FIG. 
8. 

The split process is executed by the split/merge 
processors in level 1 for backprojection. Its primary func- 
tion is to partition the projection strip into three projec- 
tion substrips according to the global location of the voxel 
25 subcubes processed by its son nodes. As mentioned herein- 
above, an overlapping region exists between projection strips 
that corresponds to two adjacent voxel subcubes. Projection 
data in the overlapping region are duplicated in two adjacent 
projection strips, as represented in FIG. 9. Before perform- 
30 ing projection splitting, the processor reads in projection 
strips from its parent node and reads in information from its 
son nodes about the starting point and dimensions of the 
voxel sector of its son nodes. Then the input projection 
strip is divided into three projection substrips, according 
35 to the global location of the voxel subcubes processed by its 
son nodes and the scanning geometry, and the resulting three 
projection strips are distributed to its son nodes. 
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The backproject process is executed by the backpro- 
ject/reproject processors in level 0 for backproject ion . Its 
primary function is to read in a projection strip for every 
view angle and backproject it to a voxel subcube. The back- 
5 projection procedure uses a voxel-driven approach and is 
illustrated as follows: 

for every view angle Q\ 

read in projection strip from its parent node 
for each voxel v in the voxel subcube 
1 0 calculate pixel location in the projection strip influenced by v for each 

pixel p influenced by v 

calculate h(v t p), the weighting factor, via interpolation 
V(v) = V(v) + h(v,p)P(v) 
N(p) = N(p) + h(p,v)h(p,v) 

15 endfor 

end for 

endfor 

Considering the parallel processing approach in 
detail with reference to the geometry of FIGS. 1, 3 and 4, an 
20 (a£ xNy xNx) voxel data set can be divided into (M) voxel sub- 
cubes, each with [((a£ / M)xN% xA#)] voxels. The starting pint 
of each voxel subcube is stated as 

for m=*0, — , (M-l) . The variable m is used to denote the 
25 voxel subcube index. For each voxel subcube, its correspond- 
ing region on the projection plane that may contribute to it 
can be identified according to the scanning geometry. There 
are M projection strips corresponding to voxel subcubes . The 
starting point of the projection strip corresponding to the 
30 m th voxel subcube is stated as 
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and each projection strip contains 

projection elements. The variable A/J denotes the number of 
detector elements along the x axis of the projection coordi- 
5 nate system. The geometric relationship between the starting 
point of each voxel subcube and the starting point of its 
corresponding projection strip is illustrated in FIG. 4. 

Each voxel in a voxel subcube handled by a backpro- 
ject/reproject processor is indexed by dv>jv>kv), where 

t-a--%(jv;-i) 

10 \M ) 

Its global coordinate can be calculated as 

V r = ($+j ¥ )xR> 

The voxel resolutions along the x, y, and z axes are speci- 
fied as (RZ,Ry,l§). The integration line connecting the x-ray 
15 source and the voxel (V,,V r ,V^) for the view angle 8i can be 
stated as 

x-psmdj _ y+ftcosft _ Z 
V x -Psm9i V ) -pcos9 i ~V z 

This integrating line will intersect the detector plane at 
the location (Xj, Yd) in the projection coordinate system show 
20 in FIG. 4. 
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atVjfCosgf + V r sine f ] 
d ~ V r Gas9,-V x ^ne,+p 

d f r cos0 j -*V z siii0 l +/} 

In order to reduce computation, the (Xd, Yd) can be calculated 
as follows: 

d B 
Y _aC 

d ~~B 

5 where 

' A = ^ + t(J^ cos 0,> sin 0, ) 
< ^=5 0 -U^sin^)+y iv (^cos^) 

and 

A = S£/2£ cos0,. +Sy /?£ sin0,. 

C 0 = s£/?* 

The values A, B, and C can be calculated by incrementing 
10 their values with a constant increment. This approach sig- 
nificantly reduces the computation complexity of calculating 
Xd and Yd. Then a bi-linear interpolating scheme is applied 

to calculate the weighting factor for all the pixels influ- 
enced by the voxel [Vx, Vy,Vz) . The computation of the pixel 
15 location {Xd, Yd) is unique and applicable to both the back 
projecting and the reprojection cases. 

The bi-level tree-type parallel-proessor architec- 
ture and the four major processing tasks (split, merge, back- 
project, and reproject) are implemented on the Meiko Comput- 
20 ing Surface (T800 transputer development system) . Each back- 
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project/reproject processor of each split/merge processor 
(FIG. 5) is mapped to a T800 transputer. 

While the illustrated embodiment employs a geometry 
wherein there is overlap only between two adjacent projection 
5 strips or subsets, it will be appreciated that more than two 
projection strips or subsets may overlap. Similarly, each 
FIG. 5 split/merge processor in level 1 has one parent node 
and three son nodes, which facilitates implementation with 
transputers which have four I/O parts. However, different 

10 numbers of son nodes may be employed. Moreover, rather than 
a three-type architecture, level 1 may comprise a single pro- 
cessor which merges all projection strips or subsets in a 
single operation, or correspondingly splits the complete pro- 
jection data set into a plurality of projection strips or 

15 subsets in one operation. 

This technique can be applied to both medical and 
industrial 3D CT imaging. Because it is based on a voxel- 
driven algebraic reconstruction technique, it is particularly 
useful for reconstructing regions of interest. The paral- 
20 lei ism. provided by the Split and Merge Reconstruction method 
can be extended to large data sets. Performance improves as 
the number of backproject/repro ject processors grows. 

While specific embodiments of the invention have 
been illustrated and described herein, it is realized that 
• 25 modifications and changes will occur to hose skilled in the 
art . It is therefore to be understood that the appended 
claims are intended to cover all such modifications and 
changes as fall within the true spirit and scope of the 
invention . 
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What is claimed is: 

1 . A parallel processing method based on the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 
5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
source scanning trajectory to obtain, at each of a plurality 
of discrete source positions 8i on the source scanning tra- 
jectory, a 2D measured cc-e beam pixel data set P ir said 
10 method comprising: 

organizing the voxel data set V as a plurality m of 
independent voxel subcubes V° through V* 1 -', each voxel subcube 
including a plurality of voxels; 

initializing the voxel data set V; and 
15 performing successive iterations to* correct the 

values of at least selected voxels of the 3D voxel data set 
based on the measured cone beam data sets P if during each 

iteration correcting the values of at least the selected 
voxels for each of the discrete source positions 6/ on the 
20 source scanning trajectory by 

reprojecting each voxel subcube V° through V m ' 1 
to calculate pixel values for a corresponding 2D 
calculated projection data subset from a group of 
calculated projection data subsets P? through P"~ l 
25 of a 2D calculated projection data set Pi including 

a plurality of pixels, the projection data subsets 
Pf through P"~ l overlapping in part 

merging the calculated projection data subsets 
P? through P*~ x to calculate pixel values of the 2D 
30 calculated projection data set P/ f the step of 

merging including adding pixel values in any over- 
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lapping regions of the calculated projection data 
subsets P? through Pr~ l , 

calculating a 2D correction projection data 
set Ei including a plurality of pixels by determin- 
ing the normalized difference between the value of 
each pixel of the measured projection data set % 

and the value of the corresponding pixel of the 
calculated projection data set Pi, 

splitting the 2D correction projection data 
set Ei into a plurality of 2D correction projection 
data subsets E- through £." , ~ l corresponding to the 
voxel subcubes V° through V™" 1 , the 2D correction 
projection data subsets E? through E?~ 1 overlapping 

in part, and the step of splitting including dupli- 
cating element values for any overlapping regions 
of the correction projection data subsets E- 
through E?~ l , and 

backprojecting each correction projection data 
subsets through £*~ ! to correct voxel values in 

the corresponding voxel subcube of the plurality of 
voxel subcubes V° through W*' 1 . 

2. A method in accordance with Claim 1, wherein: 

The step of merging is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, merging groups of calculated projection data subsets 
5 from the next lower layer to present a lesser number of cal- 
culated projection data subsets to the next higher layer, the 
lowest layer merging groups of the calculated projection data 
subsets P* through P"~ l from the step of reprojecting, and the 
highest layer producing the calculated projection data set ?/; 
10 and 
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the step of splitting is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, splitting one or more correction projection data sub- 
sets from the next higher layer to present a greater number 
15 of correction projection data subsets to the next lower 

layer, the highest layer splitting the correction data set E lf 

and the lowest layer producing the correction projection data 
subsets E° through £*~ l . 

3. A method in accordance with Claim 1, which com- 
prises employing a plurality m of backpro ject/reproject pro- 
cessors corresponding to the voxel subcubes to perform the 
steps of reprojecting and backpro jecting. 

4. A method in accordance with Claim 2, which com- 
prises employing a plurality of split/merge processors con- 
nected in a tree structure to perform the steps of merging 
and splitting. 

5. A method in accordance with Claim 1, which fur- 
ther comprises, during each iteration, for each of the dis- 
crete source positions 8i: 

in conjunction with the step of reprojecting each 
5 voxel subcube V° through V™' 1 , calculating data values for a 
corresponding 2D normalizing factor data subset from a group 
of 2D normalizing factor subsets Nf through N?~ x of a 2D nor- 
malizing factor data set N\ including a plurality of pixels, 
the 2D normalizing factor subsets Nf through N"~ x overlapping 
10 in part; 

merging the normalizing factor data subsets N? 

through N?~ l to calculate element values of the 2D normalizing 
factor data set N/; and 

during the step of calculating the 2D correction 
15 data set employing each corresponding element of the nor- 
malizing factor data set element 
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6. A method in accordance with Claim 1, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through W ml i 

5 initializing the corresponding calculated projec- 

tion data subset from the group of projection data subsets P? 

through P?~ x t and 

for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 
10 adding to the value of each pixel of the cor- 

responding calculated projection data subset from 
the group of calculated projection data subsets P° 

through P"~ x influenced by the selected voxel the 
voxel value multiplied by a weighting factor h(p,v) 
15 determined for the voxel and pixel positions based 

on geometry such that pixel values are cumulated in 
the corresponding artificial projection data sub- 
set . 

7. A method in accordance with Claim 5, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through V**' 1 : 

5 initializing the corresponding calculated projec- 

tion data subset from the group of projection data subsets Pf 
through /re- 
initializing the corresponding normalizing factor 
data subset from the group of normalizing factor data subsets 
10 Nf through ; and 

for each of the selected voxels of the 3D voxel 
data set V(v) included in the particular voxel subcube, 

adding to the value of each pixel of the cor- 
responding calculated projection data subset from 
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15 the group of calculated projection data subsets 

through P"~ l influenced by the selected voxel and 
voxel value multiplied by a weighting factor h(p,v) 
determined for the voxel and pixel positions based 
on geometry such that pixel values are cumulated in 

20 the corresponding artificial projection data sub- 

set/ and 

adding to the value of each element of the 
corresponding normalizing factor data subset from 
the group of normalizing factor data subsets N? 
25 through Nf~ l the square of the weighting factor 

h(p,v) determined for the voxel and pixel positions 
based on geometry such that elecaent values axe 
cumulated in the corresponding normalizing factor 
data subset. 

8. A method in accordance with Claim 1, wherein 
the step of backpro jecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through V 7 "**, 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data subset of the group of correction projec- 
10 tion data subsets Ef through E?~ x in a pixel posi- 

tion influenced by the selected voxel multiplied by 
a weighting factor h(p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

9. A method in accordance with Claim 6, wherein 
the step of backpro jecting is voxel drivel and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through W*' 1 , 
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5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data subset of the group of correction projec- 
10 tion data subsets E? through in a pixel posi- 

tion influenced by the selected voxel multiplied by 
the weighting factor h(p,v) determined for the 
voxel and pixel positions based on geometry such 
that corrected voxel values are cumulated. 

10. A parallel processing method based on the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 
5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
circular source scanning trajectory lying in a plane perpen- 
dicular to a rotation axis passing through the object and 
perpendicular to all planes of the detector to obtain, at 
10 each of a plurality of discrete angular positions 8/ on the 
source scanning trajectory, a 2D measured cone beam pixel 
data set said method comprising: 

organizing the voxel data set V as a plurality m of 
independent voxel subcubes V° through V**!, boundaries between 
15 adjacent voxel subcubes lying in planes perpendicular to all 
planes of the detector and parallel to each other, and each 
voxel subcube including a plurality of voxels; 

initializing the voxel data set V; and 
performing successive iterations to correct the 
20 values of at least selected voxels of the 3D voxel data set 
based on the measured cone beam data sets P ir during each 
iteration correcting the values of at least the selected vox- 
els for each of the discrete angular positions 9/ on the 
source scanning trajectory by 
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25 reprojecting each voxel subcube V° through V™- 1 

to calculate pixel values for a corresponding 2D 
calculated projection data strip from a group of 
calculated projection data strips /f through P*~ l of 
a 2D calculated projection data set Pi, the projec- 

30 tion data strips P? through P"' 1 overlapping in 

part/ 

merging the calculated projection data strips 
If through P/"~ l to calculate pixel values of the 2D 
calculated projection data set Pi, the step of 
35 merging including adding pixel values in any over- 

lapping regions of the calculated projection data 
strips P? through , 

calculating a 2D correction projection data 
set Ei by determining the normalized difference 

40 between the value of each pixel of the measured 

projection data set P g and the value of the corre- 
sponding pixel of the calculated projection data 
set Pi, 

splitting the 2D correction projection data 
45 set Ei into a plurality of 2D correction projection 

data strips E- through E?~ l corresponding to the 

voxel subcubes V° through V* 1 ^ , the 2D correction 

projection data strips E? through E?~ x overlapping 

in part, and the step of splitting including dupli- 
50 eating element values for any overlapping regions 

of the correction projection data strips E? through 
E?~ l , and 

backprojecting each correction projection data 
strip through E?~ l to correct voxel values in the 

55 corresponding voxel subcube of the plurality of 

voxel subcubes V° through V* 71 ^. 
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11. A method in accordance with Claim 10, wherein: 

the step of merging is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, merging groups of calculated projection data strips 
5 from the next lower layer to present a lesser number of cal- 
culated projection data strips to the next higher layer, the 
lowest layer merging groups of the calculated projection data 
strips P° through P"~ x from the step of repro jecting, and the 
highest layer producing the calculated projection data set P/; 
10 and 

the step of splitting is organized in a tree-type 
manner as a plurality of layers and comprises, for each 
layer, splitting one or more correction projection data- 
strips from the next higher layer to present a greater number 
15 of correction projection data strips to the next lower layer, 
the highest layer splitting the correction data set £/, and 
the lowest layer producing the correction projection data 
strips E° through E?~ l , 

12. A method in accordance with Claim 10, which 
comprises employing a plurality m of backpro ject/repro ject 
processors corresponding to the voxel subcubes to perform the 
steps of repro jecting and backpro jecting. 

13. A method in accordance with Claim 11, which 
comprises employing a plurality of split/merge processors 
connected in a tree structure to perform the steps of merging 
and splitting. 

14. A method in accordance with Claim 10, which 
further comprises, during each iteration, for each of the 
discrete source positions 9// 

in conjunction with the step of reprojecting each 
5 voxel subcube V° through V 71 ^, calculating data values for a 
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corresponding 2D normalizing factor data strip from a group 
of 2D normalizing factor strips Nf through N?~ x of a 2D nor- 
malizing factor data set Ni including a plurality of pixels, 
the 2D normalizing factor strips Nf through N?~ x overlapping 
10 in part; 

merging the normalizing factor data strips Nf 

through N?~ l to calculate element values of the 2D normalizing 
factor data set A//; and 

during the step of calculating the 2D correction 
15 data set Ei, employing each corresponding element of the nor- 
malizing factor data set 

15. A method in accordance with Claim 10 f wherein 
the step of reprojecting is voxel drivel and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through W** 1 : 

5 initializing the corresponding calculated projec- 

tion data strip from the group of projection data strips P? 

through P?~ x ; and 

for each of the selected voxels of the 3D voxel 
data set V included in the particular voxel subcube, 
10 adding to the value of each pixel of the cor- 

responding calculated projection data strip from 
the group of calculated projection data strips If 

through P"' 1 influenced by the selected voxel and 
voxel value multiplied by a weighting factor h(p,v) 
15 determined for the voxel and pixel positions based 

on geometry such that pixel values are cumulated in 
the corresponding artificial projection data strip, 

16. A method in accordance with Claim 14, wherein 
the step of reprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through V**' 1 : 
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5 initializing the corresponding calculated projec- 

tion data strip from the group of projection data strips P° 
through /re- 
initializing the corresponding normalizing factor 
data strip from the group of normalizing factor data strips 
10 N? through Nf* g and 

for each of the selected voxels of the 3D voxel 
data set V(v) included in the particular voxel subcube, 

adding to the value of each pixel of the cor- 
responding calculated projection data strip from 

15 the group of calculated projection data strips /) p 

through P"~ l influenced by the selected voxel and 
voxel value multiplied by a weighting factor h(p,v) 
determined for the voxel and pixel positions based 
on geometry such that pixel values are cumulated in 

20 the corresponding artificial projection data strip, 

and 

adding to the value of each element of the 
corresponding normalizing factor data strip from 
the group of normalizing factor data strips Nf 
25 through N?~ l the square of the weighting factor 

h(p,v) determined for the voxel and pixel positions 
based on geometry such that element values are 
cumulated in the corresponding normalizing factor 
data strip. 

17. A method in accordance with Claim 1, wherein 
the step of backprojecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through V**'*, 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube , 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
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tion data strip of the group of correction projec- 
tion data strips E? through E?~ l in a pixel position 
influenced by the selected voxel multiplied by a 
weighting factor h(p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

18. A method in accordance with Claim 15, wherein 
the step of backpro jecting is voxel driven and comprises, for 
each voxel subcube of the plurality of voxel subcubes V° 
through W 1 ' 1 , 

5 for each of the selected voxels of the 3D voxel 

data set V included in the particular voxel subcube, 

adding to the voxel value the value of each 
element of the corresponding 2D correction projec- 
tion data strip of the group of correction projec- 
10 tion data strips E? through £*" 1 in a pixel position 

influenced by the selected voxel multiplied by the 
weighting factor h(p,v) determined for the voxel 
and pixel positions based on geometry such that 
corrected voxel values are cumulated. 

19. Parallel processing apparatus implementing the 
Algebra Reconstruction Technique (ART) for iteratively recon- 
structing from cone beam projection data a 3D voxel data set 
V representing an image of an object, the cone beam projec- 

5 tion data acquired by employing a cone beam x-ray source and 
a two-dimensional array detector to scan the object along a 
source scanning trajectory to obtain, at each of a plurality 
of discrete source positions 0/ on the source scanning tra- 
jectory, a 2D measured cone beam pixel data set P t , said appa- 
10 ratus comprising: 

a data memory for storing the voxel data set V 
organized as a plurality m of independent voxel subcubes V° 
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through V****, each voxel subcube including a plurality of vox- 
els; 

15 processing elements for initializing the voxel data 

set V and performing successive iterations to correct the 
values of at least selected voxels of the 3D voxel data set 

based on the measured cone beam data sets P i9 during each 

iteration correcting the values of at least the selected vox- 
20 els for each of the discrete source positions 6/ on the source 

scanning trajectory, said processing elements including in a 
level 0 a plurality m of backpro ject/repro jectors 
corresponding to the voxel subcubes, and in a level 1 at 
least one split/merge processor, 
25 said backpro ject/repro ject processors operable 

to rep ro ject each voxel subcube V° through V™- 1 to 

calculate pixel values for a corresponding 2D cal- 
culated projection data subset from a group of cal- 
culated projection data subsets P? through P"~ x of a 
30 2D calculated projection data set Pi, the projec- 

tion data subsets ff through P"~ l overlapping in 
part, 

said at least one split/merge processor opera- 
ble to merge the calculated projection data subsets 
35 /f through P"~ x to calculate pixel values of the 2D 

calculated projection data set Pi, adding pixel 
values in any overlapping regions of the calculated 
projection data subsets P? through P*~ x , 

said processing elements including means for 
40 calculating a 2D correction projection data set Ei 

by determining the normalized difference between 
the value of each pixel of the measured projection 
data set Pi and the value of the corresponding pixel 
of the calculated projection data set Pi, 

45 said split/merge processor operable to split 

the 2D correction projection data set Ej into a 
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plurality of 2D correction projection data subsets 
Ef through E?~ l corresponding to the voxel subcubes 
V° through V* 1 " 1 , the 2D correction projection data 
50 subsets E- through E*~ l overlapping in part, dupli- 

cating element values for any overlapping regions 
of the correction projection data subsets £/ 

through E*" 1 , and 

said backproject/reproject processors operable 
55 to backproject each correction projection data sub- 

set Ef through E"" 1 to correct voxel values in the 

corresponding voxel subcube of the plurality of 
voxel subcubes V° through W 1 ' 1 . 

20, Apparatus in accordance with Claim 19, which 
comprises a plurality of split/merge processors organized as 
a tree-type structure in a plurality of layers; 

said split/merge processors operable to merge 
5 groups of calculated projection data subsets from the next 
lower layer to present a lesser number of calculated projec- 
tion data subsets to the next higher, layer, the lowest layer 
merging groups of the calculated projection data subsets ff 

through Ff 1 from the backproject/reproject processors, and 

10 the highest layer producing the calculated projection data 
set Pi; and 

said split/merge processors operable to split one 
or more correction projection data subsets from the next 
higher layer to present a greater number of correction pro- 
15 jection data subsets to the next lower layer, the highest 
layer splitting the correction data set £/, and the lowest 

layer producing the correction projection data subsets £/ 

through E*' 1 . 

21. Apparatus in accordance with Claim 20, wherein 
at least said split/merge processors comprise transputers. 
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22. Apparatus in accordance with Claim 20, wherein 
said split /merge processors of each layer each include one 
I/O port connected to the next higher layer and three I/O 
ports connected to the next lower layer. 
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